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Abstract. We analyze the thermodynamics of the atomic and (nematic) pair 
superfluids appearing in the attractive two-dimensional Bose-Hubbard model with a 
three-body hard-core constraint that has been derived as an effective model for cold 
atoms subject to strong three-body losses in optical lattices. We show that the thermal 
disintegration of the pair superfluidity is governed by the proliferation of fractional half- 
vortices leading to a Bcrczinskii-Kostcrlitz-Thouslcss transition with unusual jump 
, in the helicity modulus. In addition to the (conventional) Berczinskii-Kostcrlitz- 

Thousless transition out of the atomic supcrfluid, we furthermore identify a direct 
thermal phase transition separating the pair and the atomic supcrfluid phases, and 
^ \ show that this transition is continuous with critical scaling exponents consistent with 

those of the two-dimensional Ising universality class. We exhibit a direct connection 
between the partial loss of quasi long-range order at the Ising transition between the 
two superfluids and the parity selection in the atomic winding number fluctuations 
that distinguish the atomic from the pair supcrfluid. 
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1. Introduction 

Since the seminal work by Hohenberg [1] and by Mermin and Wagner [2], it is well known, 
that true long-range order is destroyed at finite temperatures in two-dimensional short- 
range-interacting systems with a continuous symmetry, due to of the proliferation of 
order parameter fluctuations induced by low energy excitation modes. However, quasi- 
long range order, characterized by an algebraic decay of correlations, can still persist 
at low temperatures, as is the case e.g. in the two-dimensional XY model, modelling 
the physics of e.g. coplanar magnets with U(l) symmetry. For the finite temperature 
properties of the two-dimensional XY model, vortex and anti-vortex configurations, 
constituting topological defects in the field configuration, around which the XY order 
parameter acquires a phase difference of 27rz/y with vy = ±1 (cf. the left panel of 
Figure [1]) are essential, as their unbinding beyond the Berezinskii-Kosterlitz-Thousless 
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Figure 1. The left (right) panel shows a (half-)vortex-antivortex pair in a classical XY 
model, where the phase is represented by the direction of the shown arrows. Whereas 
the arrows along the line separating the vortex-antivortex cores (denotes in the figure 
by + and -) are aligned, they form a domain wall (string) of anti-parallel spins in the 
case of half- vortices. 

(BKT) transition temperature [31 H] destroys the quasi-long-range order of the low- 
temperature phase. An even richer phase structure can be found in generalizations of 
the XY model: Directly coupling to the vortex core energy upon introducing additional 
plaquette-couplings [5] or by introducing a nematic nearest-neighbor coupling term, 
half- vortices with fractional vorticity vy = ±1/2 can be stabilzed [6j[7J|8]. Half- vortices 
exhibit a phase winding of 7r in the XY order parameter, and a half-vortex antivortex 
pair is connected by a string with a finite tension formed by a domain wall, cf. the right 
panel of Figure [TJ For sufficiently weak string tension, a BKT transition of half- vortices 
can occur upon cooling the system down from high temperatures. Below this transition 
temperature, the domain-wall strings connecting the now confined half-vortex antivortex 
pairs form. They disappear at even lower temperatures, once the strings' tension 
dominates over their entropic contribution to the free energy. This second transition 
relates to a discrete Z2 symmetry distinction between the low temperature phase and 
the (quasi-) nematic intermediate phase, and belongs to the universality class of the 
two-dimensional Ising model [61 [7J |8]. Recently, classical XY models with generalized 
nematic couplings have attracted renewed interest [HI \TU\ HTJ [12] not only because of 
fractional vortex unbinding at high temperatures but also since the transitions between 
different quasi- long-range-ordered phases show intriguing critical properties pUl [HJ EE] . 

Turning towards quantum systems, the conventional BKT-physics finds various 
realizations, e.g. by the phase fluctuations of a two-dimensional Bose gas near the 
superfluid-to-normal transition [131 E] m Helium films and in two-dimensional quantum 
gases [15], where topological vortices have been directly imaged [161 113 (HI HI 120] - 
Furthermore, half- vortex excitations have been predicted [211 E21 E31 El] and were indeed 
observed [25] in exciton-polariton condensates. Spinor condensates in the presence 
of spin-orbit coupling are also expected to give rise to the formation of half-vortex 
states [26]. Recently, attention has focused on gases of attractive bosons, where the 
above scenario of half-vortex topological defects and strings is realized due to the 
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Figure 2. Finitc-tcmpcrature phase diagram of constrained bosons on the square 
lattice for t/\U\ = 0.08 and fj,/\U\ = -0.5. The ASF-DSF transition belongs to the 
two-dimensional Ising universality class (dashed line) . The blue and red lines indicate 
the BKT transitions of the ASF and DSF, respectively and are distinguished by the 
unbinding of vortices (vy = 1) and half-vortices (yy = 1/2). The insets show different 
possible scenarios for the region where the two BKT lines and the Ising transition line 
meet. 



emergence of two different kinds of superffuid phases: besides a conventional atomic 
superffuid, a molecular superffuid of boson pairs (dimers) can form, which provides half- 
vortex topological defects [271 EE] • To prevent the attractive Bose gas from collapse, a 
scenario has been proposed recently, wherein the Bose gas is constrained to an optical 
lattice and furthermore subject to strong three-body losses that project out triple 
occupancy on each lattice site [29J. The dynamical suppression of triply occupied sites 
has been demonstrated eventually in Cesium condensates [30]. The effective model that 
describes this constrained lattice boson system is the Bose-Hubbard model with a onsite 
three-body constraint, [b\) 3 = 0, given by the Hamiltonian 

# = -t^(^+h.c.)-t'^^ 

(ij) (ij) i i 

where bi (b\) are bosonic annihilation (creation) operators and = b\bi. Here, t and 
t' denote hopping terms of atoms and dimers between nearest-neighbors (ij), U < is 
the attractive on-site interaction and the filling is controlled by a chemical potential \i. 
Furthermore, the bosonic operators have to fulfill the constraint (b\) 3 = 0, restricting 
the local occupation to two at most. For t! = 0, this model has been subject to extended 
analytical [291 ED E21 E31 [M] and numerical [351 EEl E7] studies. We are interested here 
in particular in the case of a two-dimensional square lattice geometry. Besides the 
atomic superffuid phase (ASF) with an atomic quasi-condensate (showing an algebraic 
decay of the atomic Green's function (b\bj)), a pair superffuid phase proliferates in the 
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strong coupling regime (t <C \U\) where the atomic correlations decay exponentially, but 
the boson pair correlation function {{b\) 2 (bj) 2 ) still decay algebraically. This distinct 
superfluid phase is also called the dimer superfluid phase (DSF). Recent quantum Monte- 
Carlo (QMC) calculations exhibited, that the finite temperature superfluid transition 
to the high-temperature phase is governed by the unbinding of half- vortices in the case 
of the DSF [35| [36] such that this system provides an interesting candidate for the 
investigation of such topological excitations in state-of-the-art experiments. The DSF 
phase is distinct from the ASF phase by a residual Z2 invariance, similar to the nematic 
phase of the aforementioned classical models [27J [28] . Here, we explore in particular the 
possibility of a direct thermal phase transition between the ASF and DSF region, and 
study its critical properties. 

For this purpose, we analyze the Hamiltonian of Equation (CQ) for which the 
ASF-DSF transition is driven explicitly upon varying the pair hopping amplitude t'. 
Such correlated hopping terms have been derived also for spin-1 bosons within an 
effective three-body constrained Hubbard model description [3E]. Using large scale 
QMC simulations, we focus here on exploring the thermal phase transitions in this 
model between the ASF, DSF and normal fluid (NF) phases. The transition from the 
ASF to the DSF phase can also be driven upon increasing the interaction strength \U\ 
instead of t', and one would expect the universal properties of the considered phase 
transitions to not depend on the choice of the driving parameter. However, since the 
QMC algorithm performs very efficiently in the presence of a finite pair hopping term 
t', we can access considerably larger system sizes than for the | [/(-driven transitions at 
if = [35]. 

After presenting methodological aspects in the following section, we reveal in 
Section 3 the presence of a continuous direct ASF-DSF thermal phase transition, and 
provide evidence, that this transition belongs to the two-dimensional Ising universality 
class, related to the residual Z2 symmetry distinction between the ASF and DSF phase. 
In particular, we identify relevant critical quantities and the related critical exponents, 
and show that the helicity modulus from the odd-parity subsector, with a scaling 
dimension of zero, can be used to assess the critical point rather accurately. In Section 
4, we then provide a detailed analysis of the ASF-NF and the DSF-NF thermal phase 
transitions. This allows us to clearly contrast the usual vs. unusual nature of the BKT 
transitions for the ASF and DSF phase, respectively from simulations on lattices with 
a linear size extending up to L = 200. The finite-temperature phase diagram, shown 
in Figure [21 summarizes the results of our analysis and will be further discussed in the 
following sections. 

2. Method 

The thermodynamic properties of the model presented in Equation (JTJ are here studied 
using large-scale QMC simulations; we employ a generalized directed loop algorithm in 
the stochastic series representation [391 HOI E] on a square lattice geometry with linear 
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system size L and N = L 2 lattice sites. In order to account for the dimer hopping 
term in Equation (JT]), we use two types of directed loops, where the worm can either 
carry a single boson or a boson pair (dimer) creation/annihilation operator. Typical 
system sizes simulated in this work range from L = 10 to 200 and allow for a precise 
finite-size analysis of the data, as discussed below. The simulations are performed at 
finite temperatures T, and with periodic boundary conditions taken along both spatial 
directions. 

To access the superfluid response of the system, we calculate the helicity modulus 



by measuring the winding number fluctuations W = (W x , W y ) along the x- and y- 
directions. This quantity relates directly to the superfluid stiffness, and thus allows to 
detect the onset of superfluidity at the BKT transition [121 SSI El- Furthermore, we 
define even and odd helicities, F ev e n /odd, that take into account only winding numbers 
from the corresponding parity sector, i.e. for Y even (Yodd) both W x and W y must be even 
(odd). These parity- restricted helicities allow for a discrimination of the ASF and DSF 
phases [121 [351 EE], as discussed below. Note that a factor 1/\U\ is included to make 
Y dimensionless in the above formula. The quasi-long-range order within the ASF and 
DSF phases are characterized by structure factors which correspond to the atomic and 
dimer condensate densities, C\ and C2, respectively, which are measured via the atomic 
and dimer equal-time Green's functions 



in the standard way during the off-diagonal directed-loop update [3U SH] and are defined 

as 



for a = 1,2. Previous numerical simulations exhibited a sampling problem of the dimer 
Green's function in the limit t' = 0, originating from the appearance of a fat-tailed 
distribution in the histogram of the observable for G2 inside the DSF phase [35]. The 
following simulations, however, are performed in a regime where the system is well 
within the ASF phase for t' = 0, such that these sampling problems could be avoided. 

3. Ising Transition 

The phase diagram for the model considered in Equation ([1]), shown in Figure [21 was 
obtained for fixed values of t/\U\ = 0.08 and fi/\U\ = —0.5. These values were chosen 
such that at low temperatures the system is still in the ASF phase for t' — (the ASF- 
DSF quantum phase transition for t' — is located at (t/\U\) c ~ 0.055 [35]), and thus 
exhibits quasi-long-range-order in both the atomic and dimer channel as well as both 
even and odd winding numbers. 



G 1 (i,j) = (b\b j ), G 2 (i,j) = ((btf(b j ) 2 ) 



(3) 




(4) 
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Figure 3. Atomic (Gi) and dimer (G2) equal-time Greens's function at the largest 
distance L/2 for different values of the dimer hopping strength t', at fixed T/t = 0.625. 



Upon increasing t'/t at low temperatures, the system eventually enters the DSF 
phase, due to the enhanced dimer hopping strength t'. Indeed, beyond a critical value 
of t'/t ~ 0.7, the atomic Green's function ceases to exhibit algebraic scaling, while the 
pair Green's function still remains critical. This can be seen in Figure El which shows 
both G\ and G2 at the largest distance L/2 as a function of the linear system size L 
for various values of t'/t around 0.7 at a temperature of T/t = 0.625. While the pair 
Green's function G 2 remains critical throughout this regime, the atomic Green's function 
G\ exhibits a clear suppression at large length scales for t'/t beyond 0.7, indicating the 
breakdown of its characteristic algebraic scaling in the ASF phase. The ASF-DSF 
transition has a slight temperature dependence (see Figure |2]) such that the system can 
be driven from one superffuid regime to the other one via controlling the temperature. 
As discussed above, one expects this finite temperature direct ASF-DSF transition to 
be within the two-dimensional Ising universality class. 

To assess the nature of the finite-T ASF-DSF transition, we first consider the 
behavior of various thermodynamic quantities across the transition line. In particular, 
Figure H] shows the energy E, filling n and total helicity Y for t'/t = 0.6875, all of 
which are seen to vary smoothly across the thermal transition located at {T/t) c ~ 0.73. 
Beyond this temperature, i.e. within the DSF region, the filling n shows weak but 
visible finite size effects, indicative of a different state as compared to the low-T phase. 
Based on the smooth behavior of all these quantities across the transition, we do not 
obtain any evidence for a first-order transition. We thus continue to analyze the direct 
ASF-DSF transition in terms of its critical properties. 

For this purpose, we first focus on the odd-parity helicity Y" dd, which exhibits 
distinct behavior in the ASF vs. DSF phase: The inset of Figure [5] shows Y ^ as a 
function of T for t'/t = 0.6875. This quantity indeed exhibits a pronounced change in 
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Figure 4. Energy E, rilling n and helicity Y as functions of the temperature T 
for t'/t = 0.6875, providing no indication for first-order behavior at the ASF-DSF 
transition at {T/t) c « 0.73, indicated by the dashed line. 
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Figure 5. Data collapse for the odd helicity Y ad for t'/t = 0.6875 using k = and 
v = 1. The inset shows F dd for different system sizes across the critical point at 
(T/t) c « 0.73. 



its finite-size scaling behavior at about T/t ~ 0.73: at lower T, it increases with system 
size, approaching an essentially T-independent finite value in the thermodynamic limit, 
while at larger T, it scales to zero in the thermodynamic limit. Furthermore, y o dd 
shows a crossing point of the finite-size data at criticality. This is reminiscent of the 
general scaling behavior of helicities (or superfluid densities) at second-order thermal 
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Figure 6. Data collapse for the condensate density G\ for t' /t = 0.6875 using standard 
Ising critical exponents (j3 = 1/8 and v = 1). Inset: Rcscaled condensate density C\ 
in the vicinity of the critical point exhibiting its algebraic decay C\ oc L" 1 / 4 . 



phase transitions, i.e. F dd oc L 2 ~ d . In two dimensions (d = 2) Y oc \d thus has scaling 
dimension zero, similar to the Binder ratio jlTJ 08], and can hence be used for a precise 
estimation of the critical point. This critical scaling behavior of y o dd is confirmed by 
applying a standard scaling ansatz 

~T — T 

(5) 



Y 



odd 



-k/u 



9 



Tr 



Here, k and v denote critical exponents and g is a scaling function. In order to extract 
the critical exponents as well as T c , we fit the finite-size data for Y Q ^ to the above scaling 
ansatz, approximating the scaling function g by a fourth-order polynomial. The error 
bars for the fitting parameters are estimated using a standard bootstrap analysis |49| . 
From this analysis, we obtain the values (T/\U\) C = 0.726 ± 0.006, v = 1.08 ± 0.09 and 
k = —0.01 ±0.04, consistent with the two-dimensional Ising critical exponent v = 1 [50J 
and with k = 0. Indeed, a robust data collapse of y o dd is observed in the main panel of 
Figure |5] using the Ising exponent v = 1 and n = 0. 

Further confirmation of Ising criticality at the finite-T ASF-DSF transition is 
obtained from analyzing the atomic condensate density C±. While in two dimensions 
at finite T, the expectation value of the atomic condensate density vanishes in both 
phases, it still exhibits characteristic algebraic scaling within the transition region. In 
particular, the inset of Figure |6] shows the finite-size behavior of C\ multiplied by L 1//4 , 
which makes evident characteristic scaling C± oc L -1 / 4 in the vicinity of the critical 
temperature. Using a standard scaling ansatz 

~T-T r 



d = L- 2 P' v h 



(6) 



with a scaling function h, the above scaling of C\ oc L 1 ^ at T 



T c is consistent with 



Half-Vortex Unbinding and Ising Transition in Constrained Superfluids 9 

0.08 
0.06 
^0.04 
0.02 



0.75 0.8 0.85 1 1.2 1.4 

T/t T/t 

Figure 7. Helicity Y as a function of temperature T/t for different system sizes for 
t' = (left panel) and t'/t = 0.875 (right panel) at t/\U\ = 0.08 and fi\\U\ = -0.5. The 
solid (dashed) line relates to the usual (unusual) universal helicity jump 2T/(n\U\), 
(8T/Mt/|)). 

the critical exponents /3 and v taking on their two-dimensional Ising values = 1/8 and 
v = 1 [50]. Performing an unbiased fit of the finite-size data to the above scaling ansatz 
for d, we obtain {T/\U\) C = 0.728 ± 0.009, = 0.125 ± 0.005 , and v = 0.95 ± 0.07. 
These values are indeed consistent with the Ising model values and the value of v is also 
consistent with the result from the scaling analysis of Y &&- The main panel of Figure [6] 
show the data collapse of Ci, using exact Ising exponents. The critical exponents found 
for this specific case of t'/t = 0.6875 are consistent also with values obtained along 
other cuts through the ASF-DSF phase transition line. We thus conclude, that at 
finite temperatures a direct ASF-DSF transition is present, that belongs to the two- 
dimensional Ising universality class. Furthermore, the odd-parity helicity Y Q dd, being a 
dimensionless quantity, provides a direct connection between the emergence of the Z 2 
symmetry above the ASF-DSF transition line and the loss of bare atomic contributions 
to the (odd) single-particle winding beyond the ASF region. 

4. BKT Transitions 

We next consider the phase transitions from the ASF and DSF phases into the NF 
high temperature regime. For both the ASF and DSF phases, the system undergoes a 
BKT transition at a finite temperature Tbkt- A standard way of estimating Tbkt for a 
conventional superfluid stems from the universal jump of the helicity modulus [51] and 
knowledge of the finite-size scaling behavior of the helicity at the BKT point, as derived 
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from the BKT renormalization group equations [52~| . 



(7) 



Here, Lq is of order unity and A — 1/vy is determined by the elementary vorticity vy 
at T = Tbkt- Previous studies of the thermodynamic properties at t' = showed that 
the ASF undergoes a conventional BKT transition with an unbinding of integer vortices 
[yy = 1), for which A — 1, whereas the DSF shows the unbinding of half-integer vortices 
ivy = 1/2) and the universal jump of the helicity modulus is thus increased by a factor 
of four (A = 4) [351 [36]. 

This can be seen also for the present case. In Figure [7J the helicity modulus Y 
is shown for various system sizes for two different values of t'/t, corresponding to the 
ASF {t' = 0) and DSF (t'/t = 0.875) phase, respectively. While for the ASF case, the 
data is consistent with the usual helicity jump of 2T/(tt\U\), the data for the DSF case 
exhibits a behavior that is consistent only with the larger helicity jump of 8T/(7r\U\) at 
the BKT transition. 

In order to confirm this behavior in more detail, we now employ the above scaling 
formula: Since this scaling form of the helicity modulus is valid only right at the BKT 
transition temperature, the identification of the best fit of the finite-size data to Equation 
(JZJ) indeed allows for a precise determination of the actual helicity jump, from which the 
vorticity vy can be inferred, as has been demonstrated previously, cf. e.g. Refs. [9] [TO]. 
As a measure of the goodness-of-fit, we use x 2 /DOF [53], where DOF denotes the 
number of degrees of freedom of the fit and 



is performed by summing over a range of finite system sizes Lj. Here, <Tj is the Monte 
Carlo error of Y(T,Li), and y fit denotes the fitting function in Equation (J7|). We first 
illustrate this procedure for the transition at t' — from the ASF to the NF. Figure M 
shows A and x 2 /DOF for if = 0, t/\U\ = 0.08 and fi/\U\ = -0.5. The pronounced 
minimum in x 2 /DOF in the best fit for A = 1 reproduces the standard value vy = 1 for 
a two-dimensional Bose-Hubbard model out of the ASF phase. In this case, we find that 
the position of the minimum in x 2 /DOF does not depend on the range of the system 
sizes that is included in the fitting procedure. 

For finite values of t'/t, however, we observe a mild dependence of the best fit 
value of A (i.e. the position of the minimum in x 2 /DOF) on the considered range of 
system sizes. We thus need to carefully assess the dependence of the fitting result, in 
particular the best-fit value of A, on the considered range of system sizes, since Equation 
(JTj) only holds within the scaling regime, i.e. sufficiently large values of L are available 
for the fitting procedure. In order to extract for a given ratio t'/t the best fit value of 
A = \/vy from the minimum in x 2 /DOF, we thus monitored the dependence of the 
fitting procedure on the employed range of system sizes. This is illustrated for two 
different cases, t'/t = 0.375 (ASF) and 0.75 (DSF), in Figure EB One observes that (i) 




(8) 
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Figure 8. A versus T/t for t' = and i/|J7| = 0.08 and n\\U\ = -0.5. The inset shows 
X 2 /DOF vs. A for the fit of the finite-size data using system sizes L = 10, 20, 30, .., 200 
to Equation (0. The goodness of the fit shows a clear minimum at A = 1 , indicating 
the unbinding of integer vortices. 



the minimum of x 2 /DOF shifts towards A = 1 (ASF) or 4 (DSF) as more and more 
of the smaller system sizes are omitted from the fit, and (ii) that the best-fit value of 
A converges eventually. Thus, the ASF and DSF phases are characterized by vy = 1 
and vy = 1/2, respectively and vy does not show crossover behavior as might have 
been inferred from a naive scaling analysis without monitoring the system-size-range 
convergence. Performing the above analysis for different values of t'/t, we eventually 
arrive at the phase boundaries for the ASF and DSF phases shown in Figure [21 It is 
noteworthy that Tbkt is suppressed in the vicinity of the ASF-DSF transition and gives 
rise to a superfluid re-entrance phenomenon where the system, upon increasing t'/t at 
fixed T ~ t, goes from the ASF to the NF phase and eventually enters the DSF upon 
further increasing t' . 



5. Conclusion 



In conclusion, we analyzed the finite-temperature attractive Bose-Hubbard model 
with a three-body hard-core constraint on the square lattice using large-scale QMC 
calculations. The transition to a (quasi-) nematic DSF is driven here by explicitly adding 
a nearest-neighbor pair hopping term, which allows for efficient simulations using a two- 
worm algorithm. The pair-hopping driven transition between the two quasi-ordered 
superfluids is found to be in the Ising universality class, accessible from the behavior 
of the dimensionless odd-parity helicity Y Q ^, having a scaling dimension of zero in two 
dimensions. For the case t' = 0, where the DSF proliferates in the strong-coupling 
regime, \U\ > 20t, Ng and Yang obtained evidence for a strong first-order direct ASF- 
DSF transition at finite temperatures [36]. In addition to the fact that their model 
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Figure 9. Goodness of the fit (x 2 /DOF) versus A using different fitting ranges for 
t'/t = 0.375 (left) and t'/t = 0.75 (right) at t/\U\ = 0.08 and (i\\U\ = -0.5. The 
minimum of the goodness, x 2 /DOF, is shifted towards A = 1 or 4 as smaller system 
sizes arc omitted from the fitting procedure. 



includes an additional nearest-neighbor repulsion, it is possible that the Ising transition 
will eventually be driven first-order as one approaches the strong-coupling regime in 
which the physics will be controlled by an effective Feshbach model [271 EEJ El] with a 
sizeable Feshbach detuning ~ \U\, as derived in Ref. [33] . 

An open question concerns the topology of the phase diagram (see Figure [2} in 
the region, where the two BKT transition lines and the Ising transition line meet. In 
the insets in Figure [21 we have drawn different possible scenarios, wherein the three 
lines either meet at a multicritcal point or where another (possibly Ising or first-order) 
transition line persists between one of the superfluid regions and the high temperature 
phase. Identifying the precise nature of the phase structure in this region is beyond the 
scope of this work and we leave this interesting issue for future research. 
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